//----------------------------------------------------------------------------
//
// (C) Copyrights MARIN/MSCN  2004 
//
//----------------------------------------------------------------------------
//
// $Source: $
//
// $State: Exp $
//
// $Locker:  $
//
// $Revision: 1.2 $
//
// $Date: 2003/03/25 12:52:01 $
//
// $Author: mscn $
//
// $Contractor: $
//
// $Project: $
//
// $Log: $
//
// $Description: $
//
// $End: $
//
//----------------------------------------------------------------------------

#ifndef __mst_Math_h__
#define __mst_Math_h__

#if defined(WIN32) 

#if !defined(NOMINMAX)
#define NOMINMAX
#endif

#define _USE_MATH_DEFINES

#endif

#include <math.h>
#ifndef WIN32
#include <values.h>
#else
#include <float.h>
#endif
#include <stdlib.h>

const float MST_HUGEVAL = 3.40282347e+37f;

/*! 
\def M_G 
\brief Acceleration of gravity 9.819 \f$ m/s^2 \f$
\ingroup math
*/
#define M_G 9.819f

/*!
\def MST_2DEG 
\brief Constant \f$ 180/\pi\f$
\ingroup math
*/
#define MST_2DEG 57.295779513082323

/*! 
\def M_2PI 
\brief Constant \f$ 2\pi \f$
\ingroup math
*/
#ifndef M_2PI
#define M_2PI 6.28318530717958647692
#endif

/*! 
\def M_4PI 
\brief Constant \f$ 4\pi \f$
\ingroup math
*/
#ifndef M_4PI
#define M_4PI 12.56637061435917295384
#endif

/*! 
\def M_8PI 
\brief Constant \f$ 8\pi \f$
\ingroup math
*/
#ifndef M_8PI
#define M_8PI 25.13274122871834590768
#endif

/*! 
\def M_12PI 
\brief Constant \f$ 12\pi \f$
\ingroup math
*/
#ifndef M_12PI
#define M_12PI 37.69911184307751886152
#endif

/*! 
\def M_16PI 
\brief Constant \f$ 16\pi \f$
\ingroup math
*/
#ifndef M_16PI
#define M_16PI 50.26548245743669181536
#endif

/*! 
\def M_PI 
\brief Constant \f$ \pi \f$
\ingroup math
*/
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/*! 
\def M_PI_2 
\brief Constant \f$ \frac{\pi}{2} \f$
\ingroup math
*/
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923
#endif

/*! 
\def M_PI_4 
\brief Constant \f$ \frac{\pi}{4} \f$
\ingroup math
*/
#ifndef M_PI_4
#define M_PI_4 0.78539816339744830962
#endif

/*!
\def M_PI_8
\brief Constant \f$ \frac{\pi}{8} \f$
\ingroup math
*/
#ifndef M_PI_8
#define M_PI_8 0.39269908169872415481
#endif

/*! 
\def M_1_PI 
\brief Constant \f$ \frac{1}{\pi} \f$
\ingroup math
*/
#ifndef M_1_PI
#define M_1_PI 0.31830988618379067154
#endif

/*! 
\def M_2_PI 
\brief Constant \f$ \frac{2}{\pi} \f$
\ingroup math
*/
#ifndef M_2_PI
#define M_2_PI 0.63661977236758134308
#endif

/*!
\def M_2_SQRTPI
\brief Constant \f$ \frac{2}{\sqrt{\pi}} \f$
\ingroup math
*/
#ifndef M_2_SQRTPI
#define M_2_SQRTPI 1.12837916709551257390
#endif

/*! 
\def M_SQRT2
\brief Constant \f$ \sqrt{\pi} \f$
\ingroup math
*/
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880
#endif

/*! 
\def M_SQRT1_2
\brief Constant \f$ \frac{1}{\sqrt{\pi}} \f$
\ingroup math
*/
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.70710678118654752440
#endif

#include "sys/sys.h"

#ifdef isnan
#define mst_isnan isnan
#undef isnan
#endif

#ifdef max
#define mst_max max
#undef max
#endif
#ifdef min
#define mst_min min
#undef min
#endif

#if !defined(MAXFLOAT) && defined(FLT_MAX)
#define MAXFLOAT FLT_MAX
#endif
/*! 
\def MST_MAXFLOAT
\brief Constant maximum single precision value
\ingroup math
*/
#define MST_MAXFLOAT MAXFLOAT

#if !defined(MAXDOUBLE) && defined(DBL_MAX)
#define MAXDOUBLE DBL_MAX
#endif

/*! 
\def MST_MAXDOUBLE
\brief Constant maximum double precision value
\ingroup math
*/
#define MST_MAXDOUBLE MAXDOUBLE

/*!
\def mst_clamp(v, lo, hi)
\brief Clamp value of \a v between \a lo and \a hi
\ingroup math
*/
#define mst_clamp(v, lo, hi) \
	(((v) < (lo)) ? (lo) : (v) > (hi) ? (hi) : (v))

namespace mst {

/*! 
\fn template<typename T> const T max(const T &a, const T &b)
\brief Return maximum value of \a a and \a b 
\ingroup math
*/
template<typename T> inline const T& max(const T &a, const T &b)
{
   return (a < b) ? b : a; 
}

/*! 
\fn template<typename T> const T min(const T &a, const T &b)
\brief Return minimum value of \a a and \a b 
\ingroup math
*/
template < typename T> 
inline const T& min(const T& a, const T& b)
{
   return (b < a ? b : a); 
}

/*! 
\fn template<class T> T maximum(const T &a)
\brief Return maximum value of \a a 
\ingroup math
*/
template<class T> T maximum(const T &a);

/*! 
\fn template<class T> T minimum(const T &a)
\brief Return minimum value of \a a 
\ingroup math
*/
template<class T> T minimum(const T &a);

/*! 
\fn template<class T> bool equal(const T a, const T b, const T tol)
\brief Test on equality of \a a and \a b using a tolerance \a tol
\ingroup math
*/
template<class T> bool equal(const T a, const T b, const T tol);

/*!
\fn template<class T> T factorial(T)
\brief Factorial template function
\ingroup math
*/
template<class T> T factorial(T);

/*! 
\fn template<class T> T sign(T)
\brief Sign template function 
\ingroup math
*/
template<class T> T sign(T);

/*! 
\fn template<class T> T sign(const T &a, const T &b)
\brief Sign template function 
\return The absolute value of \a a times the sign of \a b
\ingroup math
*/
template<class T> T sign(const T &a, const T &b)
{
	return (b >= 0) ? a : -a;
}

/*! 
\fn template<class T> T epsilon(T)
\brief Epsilon template function 
\ingroup math
*/
template<class T> T epsilon(T);

/*! 
\fn template<class T> T huge(T)
\brief Huge value template function 
\ingroup math
*/
template<class T> T huge(T);

/*! 
\fn template<class T> T eta(T)
\brief Eta template function 
\ingroup math
*/
template<class T> T eta(T);

/*!
\fn short factorial(short)
\brief short factorial specialisation function
\ingroup math
*/
inline short factorial(short v)
{
        if (v == 0) return 1;
        short r = 1;
        for (short i = 1; i <= v; i++) r *= i;
        return r;
}

/*!
\fn int factorial(int)
\brief int factorial specialisation function
\ingroup math
*/
inline int factorial(int v)
{
        if (v == 0) return 1;
        int r = 1;
        for (int i = 1; i <= v; i++) r *= i;
        return r;
}

/*!
\fn long factorial(long)
\brief long factorial specialisation function
\ingroup math
*/
inline long factorial(long v)
{
        if (v == 0) return 1;
        long r = 1;
        for (long i = 1; i <= v; i++) r *= i;
        return r;
}

/*!
\fn unsigned short factorial(unsigned short)
\brief unsigned short factorial specialisation function
\ingroup math
*/
inline unsigned short factorial(unsigned short v)
{
        if (v == 0) return 1;
        unsigned short r = 1;
        for (unsigned short i = 1; i <= v; i++) r *= i;
        return r;
}

/*!
\fn unsigned int factorial(unsigned int)
\brief unsigned int factorial specialisation function
\ingroup math
*/
inline unsigned int factorial(unsigned int v)
{
        if (v == 0) return 1;
        unsigned int r = 1;
        for (unsigned int i = 1; i <= v; i++) r *= i;
        return r;
}

/*!
\fn unsigned long factorial(unsigned long)
\brief unsigned long factorial specialisation function
\ingroup math
*/
inline unsigned long factorial(unsigned long v)
{
        if (v == 0) return 1;
        unsigned long r = 1;
        for (unsigned long i = 1; i <= v; i++) r *= i;
        return r;
}

/*! 
\fn short sign(short)
\brief short sign specialisation function
\ingroup math
*/
inline short sign(short v) { return (v>0) ? 1 : -1; }

/*! 
\fn int sign(int)
\brief int sign specialisation function
\ingroup math
*/
inline int sign(int v) { return (v>0) ? 1 : -1; }

/*! 
\fn long sign(long)
\brief long sign specialisation function
\ingroup math
*/
inline long sign(long v) { return (v>0) ? 1 : -1; }

/*! 
\fn float sign(float)
\brief float sign specialisation function
\ingroup math
*/
inline float sign(float v) { return (v>0.0f) ? 1.0f : -1.0f; }

/*! 
\fn double sign(double)
\brief double sign specialisation function
\ingroup math
*/
inline double sign(double v) { return (v>0.0) ? 1.0 : -1.0; }

/*! 
\fn float epsilon(float)
\brief float epsilon specialisation function (1.0e-5f) used by the Matrix4f class
\ingroup math
*/
inline float epsilon(float) { return 1.0e-5f; }

/*! 
\fn double epsilon(double)
\brief double epsilon specialisation function (1.0e-5) used by the Matrix4d class
\ingroup math
*/
inline double epsilon(double) { return 1.0e-5; }

/*! 
\fn float huge(float)
\brief float huge value specialisation function (MAXFLOAT)
\ingroup math
*/
inline float huge(float) { return MAXFLOAT; }
/*! \fn double huge(double)
\brief double huge value specialisation function (MAXDOUBLE)
\ingroup math
*/
inline double huge(double) { return MAXDOUBLE; }
/*! \fn float eta(float)
\brief float eta specialisation function (1.0e-15f)
\ingroup math
*/
inline float eta(float) { return 1.0e-15f; }
/*! \fn double eta(double)
\brief double eta specialisation function (1.0e-30)
\ingroup math
*/
inline double eta(double) { return 1.0e-30; }

/*!
\fn template<class Type> Type clamp(Type v, Type lo, Type hi)
\brief Template to clamp \a v between \a lo and \a hi
\ingroup math
*/
template<class Type> Type clamp(Type v, Type lo, Type hi)
{	Type a;

	if (lo > hi) { a = lo; lo = hi; hi = a; }
	return (((v) < (lo)) ? (lo) : (v) > (hi) ? (hi) : (v));
}

/*!
\fn template<class Type> Type round(Type a)
\ingroup math
\brief Template to calculate the rounded integer

These functions round \a a to the nearest integer, but round halfway
cases away from zero (regardless of the current rounding direction),
instead of to the nearest even integer like rint().
*/
template<class Type> Type round(Type a);

/*!
\fn template<class Type> int iround(Type a)
\ingroup math
\brief Template to calculate the rounded int

These functions round \a a to the nearest integer, but round halfway
cases away from zero (regardless of the current rounding direction),
instead of to the nearest even integer like rint().
*/
template<class Type> int iround(Type a)
{
	return (a>0) ? static_cast<int>(a + 0.5) : static_cast<int>(a - 0.5);
}

/*!
\fn template<class Type> long lround(Type a)
\ingroup math
\brief Template to calculate the rounded long

These functions round \a a to the nearest long, but round halfway
cases away from zero (regardless of the current rounding direction),
instead of to the nearest even integer like rint().
*/
template<class Type> long lround(Type a)
{
	return (a>0) ? static_cast<long>(a + 0.5) : static_cast<long>(a - 0.5);
}

/*!
\fn template<class Type> Type ceil(Type a)
\ingroup math
\brief Template to round \a a up to the nearest integer.
*/
template<class Type> Type ceil(Type a);

/*!
\fn template<class Type> Type floor(Type a)
\ingroup math
\brief Template to round \a a down to the nearest integer.
*/
template<class Type> Type floor(Type a);

/*!
\fn template<class Type> Type trunc(Type a)
\ingroup math
\brief Template to calculate the nearest integer not larger in absolute value
*/
template<class Type> Type trunc(Type a);

/*!
\fn template<class T> bool isnan(const T &a)
\ingroup math
\brief Template to check if parameter is not a number
*/
template<class T> bool isnan(const T &a);

/*!
\fn bool isnan(const float &a)
\ingroup math
\brief Single precision function to check if parameter is not a number
*/
inline bool isnan(const float &a) { return MST_ISNANF(a); }

/*!
\fn bool isnan(const double &a)
\ingroup math
\brief Double precision function to check if parameter is not a number
*/
inline bool isnan(const double &a) { return MST_ISNAND(a); }

/*!
\fn template<class T> bool finite(const T &a)
\ingroup math
\brief Template to check if parameter is a finite number
*/
template<class T> bool finite(const T &a);

/*!
\fn bool finite(const float &a)
\ingroup math
\brief Single precision function to check if parameter is a finite number
*/
inline bool finite(const float &a) { return MST_FINITEF(a); }

/*!
\fn bool finite(const double &a)
\ingroup math
\brief Double precision function to check if parameter is a finite number
*/
inline bool finite(const double &a) { return MST_FINITED(a); }

/*!
\fn template<class Type> Type exp(Type a)
\ingroup math
\brief Template to calculate the \f$ e^a \f$
*/
template<class Type> Type exp(Type a);

/*!
\fn template<class Type> Type log(Type a)
\ingroup math
\brief Template to calculate the \f$\ln{(a)}\f$
*/
template<class Type> Type log(Type a);

/*!
\fn template<class Type> Type log(Type a, Type b)
\ingroup math
\brief Template to calculate the \f$\sup{a}\log{(b)}\f$
*/
template<class Type> Type log(Type a, Type b);

/*!
\fn template<class Type> Type log10(Type a)
\ingroup math
\brief Template to calculate the \f$\sup{10}\log{(a)}\f$
*/
template<class Type> Type log10(Type a);

/*!
\fn template<class Type> Type sqrt(const Type a)
\ingroup math
\brief Template to calculate the \f$\sqrt{(a)}\f$
*/
template<class Type> Type sqrt(const Type a);

/*!
\fn template<class Type> Type sinh(Type a)
\ingroup math
\brief Template to calculate the hyperbolic sinus of \a a  which is
defined as \f$\sinh{(x)} = \frac{e^x-e^{-x}}{2}\f$
*/
template<class Type> Type sinh(Type a);

/*!
\fn template<class Type> Type cosh(Type a)
\ingroup math
\brief Template to calculate the hyperbolic cosine of \a a which is
defined as \f$\cosh{(x)} = \frac{e^x+e^{-x}}{2}\f$
*/
template<class Type> Type cosh(Type a);

/*!
\fn template<class Type> Type asinh(Type a)
\ingroup math
\brief Template to calculate the inverse hyperbolic sinus of \a a which is
defined as \f$arc \sinh{(x)} = \ln{(x + \sqrt{x^2 + 1})}\f$
*/
template<class Type> Type asinh(Type a);

/*!
\fn template<class Type> Type acosh(Type a)
\ingroup math
\brief Template to calculate the inverse hyperbolic cosine of \a a which is
defined as \f$arc \cosh{(x)} = \ln{(x + \sqrt{x^2 - 1})}\f$
*/
template<class Type> Type acosh(Type a);

/*!
\fn template<class Type> Type atanh(Type a)
\ingroup math
\brief Template to calculate the inverse hyperbolic tangent of \a a which is
defined as \f$arc \tanh{(x)} = \frac{(\ln{\frac{(1 + x)}{(1 - x)}})}{2}\f$
*/
template<class Type> Type atanh(Type a);

/*!
\fn template<class Type> Type tanh(Type a)
\ingroup math
\brief Template to calculate the hyperbolic tangent of \a a which is
defined as \f$\tanh{(x)} = \frac{e^x-e^{-x}}{e^x+e^{-x}}\f$
*/
template<class Type> Type tanh(Type a);

/*!
\fn template<class Type> Type sin(Type a)
\ingroup math
\brief Template to calculate the sinus of \a a (radians)
*/
template<class Type> Type sin(Type a);

/*!
\fn template<class Type> void sincos(Type a, Type &sinval, Type &cosval)
\ingroup math
\brief Template to calculate the sinus and cosinus of \a a (radians)
*/
template<class Type> void sincos(Type a, Type &sinval, Type &cosval);

/*!
\fn template<class Type> Type cos(Type a)
\ingroup math
\brief Template to calculate the cosinus of \a a (radians)
*/
template<class Type> Type cos(Type a);

/*!
\fn template<class Type> Type asin(Type a)
\ingroup math
\brief Template to calculate the arc sinus of \a a (radians)
*/
template<class Type> Type asin(Type a);
/*!
\fn template<class Type> Type acos(const Type &a)
\ingroup math
\brief Template to calculate the arc cosine of \a a (radians)
*/
template<class Type> Type acos(const Type &a);
/*!
\fn template<class Type> Type atan(Type a)
\ingroup math
\brief Template to calculate the arc tangens of \a a (radians)
*/
template<class Type> Type atan(Type a);
/*!
\fn template<class Type> Type tan(Type a)
\ingroup math
\brief Template to calculate the tangens of \a a (radians)
*/
template<class Type> Type tan(Type a);

/*!
\fn template<class Type> Type cosec(Type a)
\ingroup math
\brief Template to calculate the cosecant of \a a (radians) which is
defined as \f$ cosec(a) = \frac{1}{\sin{(a)}} \f$
*/
template<class Type> Type cosec(Type a) { return Type(1) / sin(a); }

/*!
\fn template<class Type> Type cosech(Type a)
\ingroup math
\brief Template to calculate the hyperbolic cosecant of \a a (radians) which is
defined as \f$ cosech(a) = \frac{1}{\sinh{(a)}} \f$
*/
template<class Type> Type cosech(Type a) { return Type(1) / sinh(a); }

/*!
\fn template<class Type> Type cot(Type a)
\ingroup math
\brief Template to calculate the cotangent of \a a (radians) which is
defined as \f$ cot(a) = \frac{1}{\tan{(a)}}\f$
*/
template<class Type> Type cot(Type a) { return Type(1) / tan(a); }

/*!
\fn template<class Type> Type coth(Type a)
\ingroup math
\brief Template to calculate the hyperbolic cotangent of \a a (radians) which is
defined as \f$ coth(a) = \frac{1}{\tanh{(a)}} \f$
*/
template<class Type> Type coth(Type a) { return Type(1) / tanh(a); }

/*!
\fn template<class Type> Type sec(Type a)
\ingroup math
\brief Template to calculate the secant of \a a (radians) which is
defined as \f$ sec(a) = \frac{1}{\cos{(a)}} \f$
*/
template<class Type> Type sec(Type a) { return Type(1) / cos(a); }

/*!
\fn template<class Type> Type sech(Type a)
\ingroup math
\brief Template to calculate the hyperbolic secant of \a a (radians) which is
defined as \f$ sech(a) = \frac{1}{\cosh{(a)}} \f$
*/
template<class Type> Type sech(Type a) { return Type(1) / cosh(a); }

/*!
\fn template<class Type> Type sqr(Type a)
\brief Template to calculate the \f$ a^2 \f$
\ingroup math
*/
template<class Type> Type sqr(Type a) { return a*a; }

/*!
\fn template<class Type> Type cubic(Type a)
\brief Template to calculate the \f$ a^3 \f$
\ingroup math
*/
template<class Type> Type cubic(Type a) { return a*a*a; }

/*!
\fn template<class Type> Type atan2(Type y, Type x)
\brief Template to calculate the arcus tangens (radians) of \a y / \a x
\ingroup math
*/
template<class Type> Type atan2(Type y, Type x);

/*!
\fn template<class Type> Type pow(Type y, Type x)
\brief Template to calculate \f$ x^y \f$
\ingroup math
*/
template<class Type> Type pow(Type a, Type b);

/*!
\fn template<class Type> Type fmod(Type a, Type b)
\brief Template to return the remainder of dividing \a a by \a b
\ingroup math
*/
template<class Type> Type fmod(Type a, Type b);

/*!
\fn template<class Type> Type abs(const Type &a)
\brief Template to return the absolute value of \a a
\ingroup math
*/
template<class Type> Type abs(const Type &a);

/*!
\fn template<class Type> bool abslt(Type a, Type b)
\brief Template to determine \f$|a| < |b|\f$
\ingroup math
*/
template<class Type> bool abslt(Type a, Type b);
/*!
\fn template<class Type> bool absle(Type a, Type b)
\brief Template to determine \f$|a| <= |b|\f$
\ingroup math
*/
template<class Type> bool absle(Type a, Type b);
/*!
\fn template<class Type> bool absgt(Type a, Type b)
\brief Template to determine \f$|a| > |b|\f$
\ingroup math
*/
template<class Type> bool absgt(Type a, Type b);
/*!
\fn template<class Type> bool absge(Type a, Type b)
\brief Template to determine \f$|a| >= |b|\f$
\ingroup math
*/
template<class Type> bool absge(Type a, Type b);
/*!
\fn template<class Type> bool nan(Type a)
\brief Template to test if \a a is not a number
\ingroup math
*/
template<class Type> bool nan(Type a);

/*!
\fn float sqrt(const float &a)
\brief float function to calculate the \f$\sqrt{a}\f$
\ingroup math
*/
inline float sqrt(const float &a) { return MST_SQRTF(a); }
/*!
\fn double sqrt(const double &a)
\brief double function to calculate the \f$\sqrt{a}\f$
\ingroup math
*/
inline double sqrt(const double &a) { return MST_SQRTD(a); }

/*!
\fn float asinh(float a)
\brief float function to calculate the \f$ arc sinh(a)\f$
\ingroup math
*/
inline float asinh(float a) { return MST_ASINHF(a); }

/*!
\fn double asinh(double a)
\brief double function to calculate the \f$ arc sinh(a)\f$
\ingroup math
*/
inline double asinh(double a) { return MST_ASINHD(a); }

/*!
\fn float acosh(float a)
\brief float function to calculate the \f$ arc cosh(a)\f$
\ingroup math
*/
inline float acosh(float a) { return MST_ACOSHF(a); }

/*!
\fn double acosh(double a)
\brief double function to calculate the \f$arc cosh{a}\f$
\ingroup math
*/
inline double acosh(double a) { return MST_ACOSHD(a); }

/*!
\fn float atanh(float a)
\brief float function to calculate the \f$ arc \tanh{(a)}\f$
\ingroup math
*/
inline float atanh(float a) { return MST_ATANHF(a); }

/*!
\fn double atanh(double a)
\brief double function to calculate the \f$ arc \tanh{(a)}\f$
\ingroup math
*/
inline double atanh(double a) { return MST_ATANHD(a); }

/*!
\fn float sinh(float a)
\brief float function to calculate the \f$\sinh{(a)}\f$
\ingroup math
*/
inline float sinh(float a) { return MST_SINHF(a); }

/*!
\fn double sinh(double a)
\brief double function to calculate the \f$\sinh{(a)}\f$
\ingroup math
*/
inline double sinh(double a) { return MST_SINHD(a); }

/*!
\fn float cosh(float a)
\brief float function to calculate the \f$\cosh{(a)}\f$
\ingroup math
*/
inline float cosh(float a) { return MST_COSHF(a); }

/*!
\fn double cosh(double a)
\brief double function to calculate the \f$\cosh{(a)}\f$
\ingroup math
*/
inline double cosh(double a) { return MST_COSHD(a); }

/*!
\fn float tanh(float a)
\brief float function to calculate the \f$\tanh{(a)}\f$
\ingroup math
*/
inline float tanh(float a) { return MST_TANHF(a); }

/*!
\fn double tanh(double a)
\brief double function to calculate the \f$\tanh{(a)}\f$
\ingroup math
*/
inline double tanh(double a) { return MST_TANHD(a); }

/*!
\fn float cos(float a)
\brief float function to calculate the \f$\cos{(a)}\f$
\ingroup math
*/
inline float cos(float a) { return MST_COSF(a); }

/*!
\fn double cos(double a)
\brief double function to calculate the \f$\cos{(a)}\f$
\ingroup math
*/
inline double cos(double a) { return MST_COSD(a); }

/*!
\fn float sin(float a)
\brief float function to calculate the \f$\sin{(a)}\f$
\ingroup math
*/
inline float sin(float a) { return MST_SINF(a); }

/*!
\fn double sin(double a)
\brief double function to calculate the \f$\sin{(a)}\f$
\ingroup math
*/
inline double sin(double a) { return MST_SIND(a); }

/*!
\fn void sincos(float a, float &sinval, float &cosval)
\brief function to calculate the float values \f$\sin{(a)}\f$ and \f$\cos{(a)}\f$
\ingroup math
*/
inline void sincos(float a, float &sinval, float &cosval) 
	{ sinval = MST_SINF(a); cosval =  MST_COSF(a); }

/*!
\fn void sincos(double a, double &sinval, double &cosval)
\brief function to calculate the double values \f$\sin{(a)}\f$ and \f$\cos{(a)}\f$
\ingroup math
*/
inline void sincos(double a, double &sinval, double &cosval) 
	{ sinval = MST_SIND(a); cosval = MST_COSD(a); }

/*!
\fn float tan(float a)
\brief float function to calculate the \f$\tan{(a)}\f$
\ingroup math
*/
inline float tan(float a) { return MST_TANF(a); }

/*!
\fn double tan(double a)
\brief float function to calculate the \f$\tan{(a)}\f$
\ingroup math
*/
inline double tan(double a) { return MST_TAND(a); }

/*!
\fn float asin(float a)
\brief float function to calculate the \f$arc \sin{(a)}\f$
\ingroup math
*/
inline float asin(float a) { return MST_ASINF(a); }

/*!
\fn float acos(const float &a)
\brief float function to calculate the \f$arc \cos{(a)}\f$
\ingroup math
*/
inline float acos(const float &a) { return MST_ACOSF(a); }

/*!
\fn float atan(float a)
\brief float function to calculate the \f$arc \tan{(a)}\f$
\ingroup math
*/
inline float atan(float a) { return MST_ATANF(a); }

/*!
\fn double asin(double a)
\brief double function to calculate the \f$arc \sin{(a)}\f$
\ingroup math
*/
inline double asin(double a) { return MST_ASIND(a); }

/*!
\fn double acos(const double &a)
\brief double function to calculate the \f$arc \cos{(a)}\f$
\ingroup math
*/
inline double acos(const double &a) { return MST_ACOSD(a); }

/*!
\fn double atan(double a)
\brief double function to calculate the \f$arc \tan{(a)}\f$
\ingroup math
*/
inline double atan(double a) { return MST_ATAND(a); }

/*!
\fn float atan2(float y, float x)
\brief float function to calculate the \f$arc tan2(y,x)\f$
\ingroup math
*/
inline float atan2(float y, float x) { return MST_ATAN2F(y,x); }

/*!
\fn double atan2(double y, double x)
\brief double function to calculate the \f$arc tan2(y,x)\f$
\ingroup math
*/
inline double atan2(double y, double x) { return MST_ATAN2D(y,x); }

/*!
\fn float pow(float a, float b)
\brief float function to calculate the \f$ a^b \f$
\ingroup math
*/
inline float pow(float a, float b) { return MST_POWF(a, b); }

/*!
\fn double pow(double a, double b)
\brief double function to calculate the \f$ a^b \f$
\ingroup math
*/
inline double pow(double a, double b) { return MST_POWD(a, b); }

/*!
\fn void srand48(long int seedval)
\brief The function srand48() initializes the random number generator.
\ingroup math
*/
inline void srand48(long int seedval) { return MST_SRAND48(seedval); }

/*!
\fn double drand48()
\brief The function drand48() returns non-negative double precision 
floating point value uniformly distributed between [0.0, 1.0>.
\ingroup math
*/
inline double drand48() { return MST_DRAND48(); }

/*!
\fn void srand(unsigned int seed)
\brief Set the seed for a new sequence of pseudo-random integers to be returned by rand.
\sa rand()
\ingroup math
 */
inline void srand(unsigned int seed) { MST_SRAND(seed); }

#ifndef RAND_MAX
#define MST_RAND_MAX mst::maximum((int&)0)
#else
#define MST_RAND_MAX RAND_MAX
#endif

/*!
\fn int rand()
\brief returns a pseudo-random integer between 0 and MST_RAND_MAX
\sa srand()
*/
inline int rand() { return MST_RAND(); }

/*!
\fn double fmod(double a, double b)
\brief Double precision function calculating the remainder of dividing \a a by \a b
\ingroup math
*/
inline double fmod(double a, double b) { return MST_FMODD(a, b); }

/*!
\fn float fmod(float a, float b)
\brief Single precision function calculating the remainder of dividing \a a by \a b
\ingroup math
*/
inline float fmod(float a, float b) { return MST_FMODF(a, b); }

/*!
\fn double log(double a)
\brief Double precision function to calculate the \f$\ln{(a)}\f$
\ingroup math
*/
inline double log(double a) { return MST_LOGD(a); }

/*!
\fn float log(float a)
\brief Single precision function to calculate the \f$\ln{(a)}\f$
\ingroup math
*/
inline float log(float a) { return MST_LOGF(a); }

/*!
\fn double exp(double a)
\brief Double precision function to calculate the \f$ e^a \f$
\ingroup math
*/
inline double exp(double a) { return MST_EXPD(a); }

/*!
\fn float exp(float a)
\brief Single precision function to calculate the \f$ e^a \f$
\ingroup math
*/
inline float exp(float a) { return MST_EXPF(a); }

/*!
\fn double log10(double a)
\brief Double precision function to calculate the \f$\sup{10}\log{(a)}\f$
\ingroup math
*/
inline double log10(double a) { return MST_LOG10D(a); }

/*!
\fn float log10(float a)
\brief Single precision function to calculate the \f$\sup{10}\log{(a)}\f$
\ingroup math
*/
inline float log10(float a) { return MST_LOG10F(a); }

/*!
\fn double log(double a, double b)
\brief Double precision function to calculate the \f$\sup{a}\log{(b)}\f$
\ingroup math
*/
inline double log(double a, double b) { return MST_LOGD(b) / MST_LOGD(a); }

/*!
\fn float log(float a, float b)
\brief Single precision function to calculate the \f$\sup{a}\log{(b)}\f$
\ingroup math
*/
inline float log(float a, float b) { return MST_LOGF(b) / MST_LOGF(a); }

/*!
\fn float ceil(float a)
\brief Single precision function round \a a up to the nearest integer
\ingroup math
*/
inline float ceil(float a) { return MST_CEILF(a); }

/*!
\fn double ceil(double a)
\brief Double precision function round \a a up to the nearest integer
\ingroup math
*/
inline double ceil(double a) { return MST_CEILD(a); }

/*!
\fn float floor(float a)
\brief Single precision function round \a a down to the nearest integer
\ingroup math
*/
inline float floor(float a) { return MST_FLOORF(a); }

/*!
\fn double floor(double a)
\brief Double precision function round \a a down to the nearest integer
\ingroup math
*/
inline double floor(double a) { return MST_FLOORD(a); }

/*!
\fn float round(float a)
\brief Single precision function to calculate the rounded value
\ingroup math
*/
inline float round(float a) { return MST_ROUNDF(a); }

/*!
\fn double round(double a)
\brief Double precision function to calculate the rounded value
\ingroup math
*/
inline double round(double a) { return MST_ROUNDD(a); }

/*!
\fn double trunc(double a)
\brief Double precision function to calculate the nearest integer not larger in absolute value
\ingroup math
*/
inline double trunc(double a) { return MST_TRUNCD(a); }

/*!
\fn float trunc(float a)
\brief Single precision function to calculate the nearest integer not larger in absolute value
\ingroup math
*/
inline float trunc(float a) { return MST_TRUNCF(a); }

/*
* abs is faster than calling fabsf
* abslt etc are faster than (abs(x1) < x2)
*/

/*!
\fn char abs(const char &a)
\brief return absolute value
\ingroup math
*/
inline char abs(const char &a) { return (a < 0) ? -(a) : (a); }

/*!
\fn short abs(const short &a)
\brief return absolute value
\ingroup math
*/
inline short abs(const short &a) { return (a < 0) ? -(a) : (a); }

/*!
\fn int abs(const int &a)
\brief return absolute value
\ingroup math
*/
inline int abs(const int &a) { return (a < 0) ? -(a) : (a); }

/*!
\fn long abs(const long &a)
\brief return absolute value
\ingroup math
*/
inline long abs(const long &a) { return (a < 0) ? -(a) : (a); }

/*!
\fn float abs(const float &a)
\brief return absolute value
\ingroup math
*/
inline float abs(const float &a) { return (a < 0.0f) ? -(a) : (a); }

/*!
\fn double abs(const double &a)
\brief return absolute value
\ingroup math
*/
inline double abs(const double &a) { return (a < 0.0) ? -(a) : (a); }

/*!
\fn bool abslt(float a, float b)
\brief return \f$ |a| < b\f$
\ingroup math
*/
inline bool abslt(float a, float b) 
	{ return (a) < (b) && -(a) < (b); }

/*!
\fn bool abslt(double a, double b)
\brief return \f$ |a| < b\f$
\ingroup math
*/
inline bool abslt(double a, double b) 
	{ return (a) < (b) && -(a) < (b); }

/*!
\fn bool absgt(float a, float b)
\brief return \f$ |a| > b\f$
\ingroup math
*/
inline bool absgt(float a, float b)
	{ return (a) > (b) || -(a) > (b); }

/*!
\fn bool absgt(double a, double b)
\brief return \f$ |a| > b\f$
\ingroup math
*/
inline bool absgt(double a, double b)
	{ return (a) > (b) || -(a) > (b); }

/*!
\fn bool absle(float a, float b)
\brief return \f$ |a| <= b\f$
\ingroup math
*/
inline bool absle(float a, float b) 
	{ return (a) <= (b) && -(a) <= (b); }

/*!
\fn bool absle(double a, double b)
\brief return \f$ |a| <= b\f$
\ingroup math
*/
inline bool absle(double a, double b) 
	{ return (a) <= (b) && -(a) <= (b); }

/*!
\fn bool absge(float a, float b)
\brief return \f$ |a| >= b\f$
\ingroup math
*/
inline bool absge(float a, float b) 
	{ return (a) >= (b) || -(a) >= (b); }

/*!
\fn bool absge(double a, double b)
\brief return \f$ |a| >= b\f$
\ingroup math
*/
inline bool absge(double a, double b) 
	{ return (a) >= (b) || -(a) >= (b); }

/*!
\fn bool nan(float a)
\brief Single precision function to test if \a a is not a number
\ingroup math
*/
inline bool nan(float a) { return MST_ISNANF(a); }

/*!
\fn bool nan(double a)
\brief Double precision function to test if \a a is not a number
\ingroup math
*/
inline bool nan(double a) { return MST_ISNAND(a); }

/*! 
\fn bool equal(const double a, const double b, const double tol)
\brief double equality check with tolerance
\ingroup math
*/
inline bool equal(const double a, const double b, const double tol)
{
	return (a >= b - tol) && (a <= b + tol);
}

/*! 
\fn  bool equal(const float a, const float b, const float tol)
\brief float equality check with tolerance
\ingroup math
*/
inline bool equal(const float a, const float b, const float tol)
{
	return (a >= b - tol) && (a <= b + tol);
}

/*!
\fn float maximum(const float &)
\brief return maximum single precision floating point value
\ingroup math
*/
inline float maximum(const float &) { return MAXFLOAT; }
/*!
\fn float minimum(const float &)
\brief return minimum single precision floating point value
\ingroup math
*/
inline float minimum(const float &) { return -MAXFLOAT; }

/*!
\fn double maximum(const double &)
\brief return maximum double precision floating point value
\ingroup math
*/
inline double maximum(const double &) { return MAXDOUBLE; }

/*!
\fn double minimum(const double &)
\brief return minimum double precision floating point value
\ingroup math
*/
inline double minimum(const double &) { return -MAXDOUBLE; }

/*!
\fn char maximum(const char &)
\brief return maximum char value
\ingroup math
*/
inline char maximum(const char &) { return char(~(1 << ((8*sizeof(char))-1))); }

/*!
\fn int maximum(const int &)
\brief return maximum int value
\ingroup math
*/
inline int maximum(const int &) { return ~(1 << ((8*sizeof(int))-1)); }

/*!
\fn char minimum(const char &)
\brief return minimum char value
\ingroup math
*/
inline char minimum(const char &) { return -1 << (sizeof(char)*8-1); }

/*!
\fn int minimum(const int &)
\brief return minimum integer value
\ingroup math
*/
inline int minimum(const int &) { return -1 << (sizeof(int)*8-1); }

/*!
\fn long maximum(const long &)
\brief return maximum long value
\ingroup math
*/
inline long maximum(const long &a) 
	{ int s = 8*sizeof(a) - 1; return ~(1 << s); }

/*!
\fn long minimum(const long &)
\brief return minimum long value
\ingroup math
*/
inline long minimum(const long &a) 
	{ int s = 8*sizeof(a) - 1; return -1 << s; }

/*!
\fn short maximum(const short &)
\brief return maximum short value
\ingroup math
*/
inline short maximum(const short &) { return short(~(1 << ((8*sizeof(short))-1))); }

/*!
\fn short minimum(const short &)
\brief return minimum short value
\ingroup math
*/
inline short minimum(const short &) { return short(-1 << (sizeof(short)*8-1)); }

/*!
\fn unsigned char maximum(const unsigned char &)
\brief return maximum unsigned char value
\ingroup math
*/
inline unsigned char maximum(const unsigned char &) { return ~((unsigned char)0); }

/*!
\fn unsigned int maximum(const unsigned int &)
\brief return maximum unsigned int value
\ingroup math
*/
inline unsigned int maximum(const unsigned int &) { return ~((unsigned int)0); }

/*!
\fn unsigned char minimum(const unsigned char &)
\brief return minimum unsigned char value
\ingroup math
*/
inline unsigned char minimum(const unsigned char &) { return 0; }

/*!
\fn unsigned int minimum(const unsigned int &)
\brief return minimum unsigned int value
\ingroup math
*/
inline unsigned int minimum(const unsigned int &) { return 0; }

/*!
\fn unsigned long maximum(const unsigned long &)
\brief return maximum unsigned long value
\ingroup math
*/
inline unsigned long maximum(const unsigned long &) { return ~((unsigned long)0); }

/*!
\fn unsigned long minimum(const unsigned long &)
\brief return minimum unsigned long value
\ingroup math
*/
inline unsigned long minimum(const unsigned long &) { return 0; }

/*!
\fn unsigned short maximum(const unsigned short &)
\brief return maximum unsigned short value
\ingroup math
*/
inline unsigned short maximum(const unsigned short &) { return ~((unsigned short)0); }

/*!
\fn unsigned short minimum(const unsigned short &)
\brief return minimum unsigned short value
\ingroup math
*/
inline unsigned short minimum(const unsigned short &) { return 0; }

/*!
\enum fpclass_t
\brief Floating point value classification enumeration
\ingroup math
*/
enum fpclass_t {
	eNan,		//!< Value is not a number
	eNInfinite,	//!< Value is negative infinity
	ePInfinite,	//!< Value is positive infinity
	eNZero,		//!< Value is negative zero
	ePZero,		//!< Value is positive zero
	eNSubnormal,	//!< Value is negative too small to be represented in normalized format
	ePSubnormal,	//!< Value is positive too small to be represented in normalized format
	eNNormal,	//!< Value is a negative valid floating point value
	ePNormal	//!< Value is a positive valid floating point value
};

/*!
\fn template<class T> fpclass_t fpclass(const T &a)
\ingroup math
\brief Template to classify the floating point value
*/
template<class T> fpclass_t fpclass(const T &a);

/*!
\fn fpclass_t fpclass(const float &a)
\ingroup math
\brief Single precision function to classify the floating point value
*/
inline fpclass_t fpclass(const float &a) { MST_FPCLASS(a); }

/*!
\fn fpclass_t fpclass(const double &a)
\ingroup math
\brief Double precision function to classify the floating point value
*/
inline fpclass_t fpclass(const double &a) { MST_FPCLASS(a); }

/*!
\fn template<class T> bool infinity(const T &a)
\ingroup mat
\brief Determine if the argument is positive or negative infinity
*/
template<class T> bool infinity(const T &a)
{
	return (fpclass(a) == eNInfinite) || (fpclass(a) == ePInfinite);
}

}

#endif

